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Abstract 


A planar simulation of film boiling and bubble formation m water at 373°C, 219 bar on 
an isothermal horizontal surface was performed by using a volume of fluid (VOF) based 
interface tracking method The complete Navier-Stokes equations and thermal energy 
equations were solved m conjunction with a interface mass transfer model The numerical 
method takes into account the effect of temperature on the transportive thermal prop- 
erties (thermal conductivity and specific heat) of vapor, effects of surface tension, the 
interface mass transfer and the corresponding latent heat The computations provided 
a good insight into film boiling yielding quantitative information on unsteady periodic 
bubble release patterns and on the spatially and temporally varying film thickness The 
computations also predicted the transport coefficients on the horizontal surface, which 
were greatly influenced by the variations m fluid properties, compared to calculations 
with constant properties 
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Cp,c v specific heat at constant pressure/ volume 
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Nomenclature 


v y component of velocity 

v fluid velocity 

x spatial coordinate m the horizontal direction 

x t ,Xj space vectors 
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Chapter 1 


Introduction 


Boiling flows are ubiquitous m the energy and processing industries due to the fact that 
phase change processes are an efficient way of transporting energy Many aspects of boil- 
ing flows that are not well understood even today. The small spatial scales and fast time 
constants of many of the physical processes associated with the boiling are the imped- 
iments of reliable experimental data Advanced computational methods are capable of 
providing solutions for fluid flow problems with moving interfaces separating gas and liq- 
uid phases 


The phase change problems are dependent on the simultaneous coupling of many effects 
most of them can not be ignored The modeling of mass, momentum and energy trans- 
port must include surface tension, latent heat, interphase mass transfer, discontinuous 
material properties and complicated liquid- vapor interface dynamics Here in this work 
we have included all these properties while predicting and analyismg the formation of 
vapor bubbles m film boiling region 
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Introduction 


1.1 Description of the Problem 


Film boiling is characterized by the existence of a continuous vapor film between a heated 
surface and the liquid layer In this work, the bubble formation and heat transfer on a 
horizontal surface have been numerically analyzed using a volume of fluid (VOF) based 
interface tracking method For this, we consider a heated wall with a uniform temperature 
of 388° C, superheated by 15° C relative to the pool of saturated water at near-critical 
conditions of 373° C, 219 bar (T c = 374 15°C, P c — 221 29 bar). The model assumes that 
location of the bubbles are spaced on a solid surface m a square pattern separated by the 
Taylor fastest-growing wavelength given by Berenson [2] 


A “ =2,r \£S (ii) 

The bubble diameter and the bubble height m Berenson’s model are considered to be 
proportional to the bubble spacing. However, the characteristic length m the present 
investigation, for presenting the pertinent dimensionless parameters, such as, Nusselt 
number has been considered as 



(Pi ~ Pg) 9 


(1 2 ) 


1.2 Scope of the Present Work 


In this investigation we explore by the numerical simulations discussed in the subsequent 
chapters, the nature of the heat transfer processes and particularly the influence of the 
large variations m the thermal and transport properties of the vapor 


The present study has significant applications in order to understand thermal issues in 
high speed civil transport, heat treatment of metals, nuclear reactor design, cryogenics, 
cooling for superconducting applications and many other thermal storage systems 
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1.3 Layout of the Thesis 

In chapter- 1, we have already discussed the genesis of the problem Chapter-2 provides a 
review of the literature m the boiling flows In this chapter, different methods for simulat- 
ing free surface flows and the evaluation of volume of fluid methods (VOF) have also been 
discussed The mathematical formulation of the problem for simulation is presented m 
Chapter-3 This chapter contains the governing equations, boundary conditions, interface 
tracking algorthim, jump conditions at liquid vapor interface and numerical procedure 
Kernel function calculation for smoothing void field and LVIRA algorithm for defining in- 
terface location are also discussed here in detail Chapter-4 present results for the current 
simulation Here the nature of periodic bubble release pattern, distribution of wall heat 
flux at various time instant, temperature contours and velocity vectors are shown and dis- 
cussed The comparison of computed space averaged Nusselt number value with different 
correlations is also done m this chapter. Chapter-5 includes the concluding remarks and 
the scope for the further research 



Chapter 2 


Review of Literature 


2.1 Review of Boiling Flows 


Numerous publications m the open literature document that there has been a long-lasting 
interest among the researchers to clarify interface transport mechanisms m liquid-vapor 
phase change processes In the past, appropriate studies m this field were carried out 
experimentally Albeit having very useful contributions m the development of the sub- 
ject, the early investigations could not provide the physical details needed for a closer 
understanding of the bubble formation and the time varying heat transfer characteristics 
Experimental studies m boiling have yielded several empirical correlations that are valid 
for specific cases Only the advancement of computational techniques have opened up 
new ways to carry out the investigations m detail on boiling, revealing the nuances of the 
interface transport processes. The volume of fluid (VOF) method of Hirt and Nichols [9] 
forms the building block of computations involving two fluids separated by a sharp in- 
terface The VOF method has been modified very successfully by Welch and Wilson [30] 
to simulate horizontal film boiling Son and Dhir [22] have also performed complete nu- 
merical simulation of the evolution of the vapor-liquid interface during film boiling on a 
horizontal surface. In another approach, Son and Dhir [23] have modified the level set 
method of Sussman et al [24] to perform numerical simulation of film boiling At low wall 
superheats, they have observed the upward movement of the interface initiating the vapor 
bubble formation Having released the bubble, the interface drops down alternatively at 
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the nodes and anti-nodes of a Taylor wave 


June and Tryggvason [11] have performed excellent simulations of film boiling using single 
field formulation where one set of conservation equations are written for the entire flow 
field and different phases are treated as one fluid with variable material properties They 
have used source terms m the continuity and the energy equations as an enhancement of 
the method of Unverdi and Tryggvason [25] Interfacial source terms for surface tension, 
interface mass transfer and latent heat are added as delta functions that are non-zero 
at the interface boundary Banerjee and Dhir [1] have performed direct simulation of 
evolving interface during sub-cooled film boiling The simulation provides the shape and 
growth rates of the interface and the associated thermal behavior 


The overview on numerical simulations of pool boiling has been authoritatively reviewed 
by Dhir [7] He pointed out that numerical simulations of evolving liquid-vapor interfaces 
during a phase change process such as boiling provide significant additional insights into 
the phenomena More recently Welch and Rachidi [28] developed another VOF based 
interface tracking method to explore film boiling on a horizontal surface and considered a 
conjugate problem between a solid wall and the boiling fluid They carried out the analy- 
sis considering constant thermal properties at the saturation temperature, whereas it has 
been seen from the NBS (National Bureau of Standards /NRC (National Research Coun- 
cil of Canada) Table (Table 4 1) that at near critical pressure, the transportive thermal 
properties of vapor vary significantly with temperature The present work is an exten- 
sion of the work by Welch and Rachidi [28] with respect to the variable thermal properties 


2.2 Overview of Methods for Simulation of Free 
Surface Flows 


Fluid flows with free surfaces or material interfaces occur in a large number of natural 
and technological processes Free surfaces mean here to be the surfaces on which dis- 
continuities exist m one or more variables Often properties that depend on the shape 
of the interface itself (e g , surface tension) play an important role m the dynamics of 
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the problem, and the physics (e g , capillarity) can drive the flow making it essential to 
accuiately detennine 'he position, cuivatuie, and topology of the mteiface 

All numerical simulations of fiee suiface flows lequne a discietisation 01 nodalization 
to allow numerical treatment on computeis There are two fundamentally different ap- 
pioaches Eulenan methods use a frame of reference (discretisation gild 01 mesh, contiol 
volumes, nodes, etc ) fixed in space, and matter (fluids) moving through this frame of 
reference Lagrangian methods instead use a frame of reference (marker paiticals, grid 
01 mesh, etc ) fixed to and moving with the matter Lagrangian method is less useful if 
strong distortions of the flow occur, or if the topology changes, foi example, as rno\ mg 
boundaries merge or break up Then, the frame of reference has to be re-imtialized (1 en- 
meshing, re-ordermg of points, etc ), which may be very cumbersome Eulenan methods, 
on the other hand, tend to be less accurate due to numerical diffusion 

For our problem, we expected strong distortions, and many topology changes during 
detachment or break-up of bubble We therefore chose an Eulerran method, m spite of its 
lower accuracy 

The first method capable of modeling gas-liquid flow, separated by a moving interface, 
was the Markei and Cell method (MAC) of Harlow and Welch [ 8 ] This was m fact 
a combination of an Eulenan solution of the basic flow field, with Lagiangian maiker 
particles attached- to the liquid to distinguish it from the gas While the staggered mesh 
layout and other features of MAC have become a model for many othei Eulenan coades, 
the marker particles proved to be computationally too expensive and have been raiely 
used. 

Subsequent to MAC, Volume of Fluid (VOF) method of Hirt and Nichols [ 9 ] becane 
popular Instead of marker particles, a liquid volume-fraction field, usually named a is 
used It indicates for each computational cell its liquid fraction (a = 1 liquid, a = 0 
gas, and 0 < a < 1 for an interfacial cell), a is a property of the fluid with which it 
moves The key feature of a VOF method lies m its technique for advecting a standard 
finite-difference methods suffer from numerical diffusion, which smoothen a excessively 
Therefore, special, often geometry-based, methods are used to advect a 

The a field is the only phase information stored m VOF methods Approximate interface 
locations are found from a so-called mteiface reconstruction This is needed foi advect- 
ing a, for determining the local properties (density, viscosity) and for better graphical 
representation In earlier versions, usually called SLIC (Simple Line Interface Method) 
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after Noh and Woodward [15], the interface is approximated as piecewise constant, where 
interface within each cell are assumed to be the lines aligned with one of the logical mesh 
coordinates 

More accurate reconstructions are possible with PLIC (Piecewise Linear Interface Con- 
struction) methods, pioneered by Youngs [31] Youngs’ method positioned each recon- 
structed interface line, defined by a slope and intercept, within the control volume cell 
The slope of the line is given by the interface normal (gradient of the volume fractions) , 
and the intercept follows from invoking the volume conservation Many extensions and 
enhancements to the significant work of Youngs have occurred since its introduction Re- 
cently, modification of Youngs method has been done with the help of Least Squares Vol- 
ume of Fluid Interface Reconstruction algorithm (LVIRA) method of Puckett et al [17] 
This method is well implemented by Welch and Rachidi [28] for exploring film boiling on 
a horizontal heated surface. This method has several advantages the fluxes of a , with 
which the phase field a is updated, can be determined more accurately, and essentially free 
of numerical diffusion The interface area is known approximately, so that surface tension 
and mass transfer can be applied more easily Fluid properties (density and viscosity) can 
be allocated accurately Finally, the straight lines also give a graphical representation of 
good quality 

Most earlier versions of PLIC used rectangular, regular grids, and were in two dimensions 
only Kothe et al. [13] and Welch et al. [29] showed that PLIC-VOF can be extended to 
general, hexahedral meshes in three dimensions without additional difficulty 

The Level Set Method of Osher and Sethian [16] is another useful method for interface 
tracking and reconstruction For interface tracking, they use a new field variable (j) (a 
smooth function), the absolute value of which indicates the shortest distance to the in- 
terface, but is positive in one fluid and negative in the other The interface is at the set 
of points where <f> = 0 

Much the same as a in the VOF method, 4> is approximated as a property of the fluid 
with which it moves, and thus follows a simple advection equation Dtfi/Dt — 0 However, 
due to flow distortions, <j> deviates more and more from the true distance to the interface, 
and therefore has to be re-initialized every couple of time-steps by the solution of a 
differential equation. For numerical reasons, properties such as density or viscosity cannot 
change sharply at the interface, but require a smooth transition (usually fitted with a sine 
function) of about three to five meshes wide Although, 0 can be discretely conserved, 
the mass enclosed by the zero level set of (j) is not conserved However, Chang et al [5] 
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True Interface 


sue 


PLIC 



Figure 2 1 A true interface, and the corresponding void field. Reconstruction of the 
interface by horizontal or vertical lines (SLIC) and by straight lines of arbitrary inclination 
(PLIC) 
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that this can be overcome 


2.3 Epilogue 


The Volume of Fluid (VOF) and the Level Set Method are the two attractive methods 
for handling flows with interface and free surface However, we preferred to use a variant 
of Volume of Fluid method The method followed here is computationally simple for its 
explicit way of advancing the liquid- vapor interface Furthermore, the method includes a 
very elegant version of a Continuum Surface Tension (CST) model of Brackbill et al , [3] 
The CST model will be discussed m a subsequent chapter The numerical technique 
m this investigation also uses a modified version of LVIRA method due to Welch and 
Rachidi [28] 



Chapter 3 


Formulation of the Problem 


3.1 Governing Equations 


The mass, momentum and energy conservation equations for the incompressible Newto- 
nian fluids for the liquid and vapor phases are given by 


(au, au t v,\ dp a ( JU,\ 

\ dt dx t ) dxj dx l \ dx l ) 

( d (P c P g ) + d(p Cp9 Uj) \ _ _d_f k dj\ 
V dt dxj J dxj y dxj J 


( 31 ) 
(3 2 ) 
(3 3 ) 


Here, U 3 , p, Cp, p, 9, p, and k are the fluid velocity, pressure, specific heat, density, 
temperature, viscosity and thermal conductivity, respectively The dissipation term in 
the energy equation has been neglected 
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Figure 3 1 The domain for the simulation of film boiling 


3.2 Boundary Conditions 


Figure 3 1 shows the domain of interest for the present investigation The simulation is 
planar, two-dimensional The computational domain is given by ABCD The boundary 
conditions are symmetry conditions at the left and right boundaries 


dv 00 

at x = 0 and x — A 0 /2 « = 0, — = 0, — = 0 

ox ox 


Constant wall temperature condition is used on the solid-fluid interface 


at y — 0 0 — 0 sup 
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Outflow boundary conditions are used on the top surface of the domain 


at y — H 


n p-p 

o — o — a — — Jr c 

oy ay dy 


The outlet pressure is the saturation pressure less the hydrostatic pressure difference from 
the initial film level to the outlet The boundary condition at the vapor liquid interface is 
of special concern m this study In order to address this issue, a suitable interface tracking 
method has been implemented 


3.3 Interface Tracking 


The presence of two phases m the fluid requires handling of the phase interface We 
advect the interface using enhancement of VOF method of Hirt and Nichols [9] due to 
Youngs [31] The method of Youngs is implemented at the end of a time cycle to calculate 
the new density field using conservation of mass for each cell 

r\ 

— f pdV + / pv ■ ndS = 0 (3 4) 

at JVc JSc 

where V c is the cell volume and S c is the cell surface The symbol v is used for the fluid 
velocity Once, the new cell densities are found, the cell void fractions are calculated using 


a = 


P-P<? 
Pi ~ Pg 


(3 5) 


Here pi and p g are the densities of the saturated liquid and saturated vapor, respectively 
The implementation of the method of Youngs has been well documented m Rider and 
Kothe [20], Rudman [21] and Welch and Wilson [30] The phase interface is represented 
by piecewise linear curves Figure 3 2 shows a typical two phase cell The orientation of 
the curve within each two phase cell is determined using the procedure based on LVIRA 
method of Puckett et al [17] 
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3.3.1 Least Squares Interface Reconstruction Algorithm 
(LVIRA) 


Least Squares Interface Reconstruction Algorithm (LVIRA) is a volume-of-fluid interface 
reconstruction algorithm for constructing an approximation to the interface at the be- 
ginning of every time cycle, given the values of the void fraction a tJ This algorithm 
produces a linear approximation to the interface which are also continuous across the cell 
boundaries, m each two phase cell The interface is uniquely defined by the unit normal, 
n, the orientation of the interface and the offset distance, l Therefore, LVIRA can be 
understood as a function which returns n, l, given the void fraction of each two-phase cell 
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In the LVIRA the approximate interface in a two phase cell ( 1 , j) is determined by mini- 
mization of the following function in the 3x3 block of cells centered on the ( 1 , j) cell 

i 

Gv?( n ) = X] {&i+k,j+i ~ d(n, (3 6) 

subject to the constraint at the center cell of the three by three block of cells 


& ( n > 0* J 


(3.7) 


Here is the actual void fraction of the cell (i, j) and &(n ,l) h0 is the function that maps 
the line with orientation n and offset length l into a void fraction of the cell (i,j) The offset 
length, l, is obtained by satisfying the constraint The function G v (n) is the squared 
deviation of the actual void fraction from the void fraction given by the mapping using the 
same line for the entire three by three block of cells centered at cell (i,j). Minimization of 
this function rotates and translates the line m such a way as to ensure that the mapping 
is exact for the center cell and that the straight line associated with this mapping is the 
best fit to the void field for the neighboring cells. The result of the minimization is the 
calculation of the unit normal and the offset length ( n, l) for the cell at the center of the 
three by three block This calculation is performed for all the mixture cells, satisfying the 
constraints that the mapping d(n, I) tJ returns the actual void fraction This constr amed 
minimization of the function G X] (n) is a nonlinear problem requiring an initial estimate 
for n This initial n is the value determined by the modified VOF method of Youngs [31] 
as 


Vo: 

jv^i 


(3 8) 


We note that this algorithm is second order method and has the property that it 
reproduces linear interfaces exactly 


3.3.2 Direction Split Approach 

Given the orientation of the planar surface that represents the interface in a cell, we 
determine the location of the oriented surface, such that the surface partitions the cell 
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Figure 3 3 Schematic of cell flux calculation 
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into liquid and gas regions of correct volume, based upon the void fraction (which is defined 
as the volume fraction of the liquid) of the cell Given the location of the planar interface 
(1 e , knowing n, l) in each two phase cell and the velocity at a junction face between two 
cells, the mass flux is determined from simple volumetric considerations Once the mass is 
fluxed across the cell boundaries m one direction (Fig 3 3), the interface is reconstructed 
before mass is fluxed m the second direction This is the Direction Split Approach of 
Rudman [21] 


3.4 Jump Conditions at Liquid Vapor Interface 


The mass transfer across the interface is modeled following Welch and Rachidi [28] We 
consider a computational cell containing a volume of the liquid phase adjacent to a volume 
of the vapor phase We write a mass balance for each phase as 


-7- [ pdV + [ pv ndS + f p(v — vA ndS = 0 

dtJvgit) Js g (t ) 

-7- [ pdV + [ pv ndS + [ p(v — vj) ndS = 0 

dt JViit) Js,(t) J Sj(t) 


(3 9) 
(310) 


Here, VJ, Si, V g , S g are the volume and surface of the liquid and vapor regions, respectively 
Sj is the phase interface at the common boundary of the two regions, moving with velocity 
v/ The normal vector n points into the liquid phase on 5; From the above, and 
taking into account the incompressibility of each phase, and under the situation that the 
overall volume is time invariant, the conservation of mass statement for the cell volume 
is determined ' 



(311) 


Here, ||.|| indicates the jump in the variable of interest across the phase interface and S c 
is the surface bounding the computational cell The mass and energy jump conditions at 
the interface may be estimated as 


(3 12) 


||p(v-v,)|| n = 0 
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\\ph{v- V/ )|| *n = -||q|| n (3 13) 

entails the contribution of jump m the conservation of mass equation as 

||(v-vl)|| n=(i-L)M£ (3 14) 

Pi Pg rllg 

Here, h is the enthalpy and h[ g = h g — hi is the latent heat of vaporization while q is the 
heat flux vector We assume the phase interface to be at the saturation temperature of 
the liquid pressure 

Oi = 9 S at(Pi) (3 15) 

The set of equations, Eq (3 13) through Eq (3 15), incorporates various simplifying 
assumptions Surface properties other than surface tension are neglected We neglect 
kinetic energy and viscous work terms as well as surface tension work terms in the energy 
jump and we have also neglected the viscous dissipation m the energy equation These 
are common approximations in the analysis of liquid- vapor phase change phenomena [4] 
and have been used m the numerical studies of liquid-vapor phase change [22, 23, 27] 
The temperature condition given by Eq (3 15) is a common approximation for which 
justification may be found m [23] 

The energy jump condition includes the effect of heat flux vector (normal to the inter- 
face) on the mass transfer rate across the interface The heat flux vector will generally be 
discontinuous and any smoothing of this vector will distort the mass transfer amount We 
utilize the interface geometry associated with the VOF method to construct a proper heat 
flux jump source term for use m Eq (3 14) By a proper heat flux jump, we mean that the 
normal components of the temperature gradients are calculated without reaching across 
the phase interface The required geometry is provided m the advection step described 
above and is shown in Fig 3 3 Given the unit normal vector, n and the parameter l which 
provides the location of the interface we apply the temperature condition, Eq. (3 15), to 
the point located at the center of the piecewise linear segment (the interface point m Fig 
3 2) It is then a simple matter to calculate liquid side temperature gradients as well as 
vapor side temperature gradients We then are able to construct a proper heat flux jump 
across the phase interface by multiplying the normal temperature gradients in each phase 
by their corresponding conductivities In addition, it is important that the cells proximal 
to the interface but not containing the phase interface see proper temperature gradients 
m the convection and diffusion terms of the energy equations We ensure this by using the 
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temperature gradients used m forming the energy jump condition to extrapolate liquid 
and vapor temperatures at mixture cell centers thus ensuring that the cells neighboring 
mixture cells also see the proper temperature gradients 


3.5 Modified Momentum Equation 


The momentum equations are augmented using the continuum surface tension model of 
Brackbill et al [3] m the following way 

Q\ r rp 

/)(«)(— + v Vv) = -VP + />(d) g + V [/i(a)(Vv + (Vv) J )] + aaVa (3 16) 

where a is a smoothened void field and k is the curvature of the surface defined by 
smoothened void field The surface tension force is applied to a transition region of nine 
cells thick centered at the interface The density and viscosity vary with the void fields 
as 

p(a) — apt + (1 — a)pg (3 17) 

p(a) = ani + (1 - a) p.g (3 18) 

In the sequel, it can be said that for the interface cells, we use the augmented momentum 
Eq (3 16), the modified conservation of mass Eq. (3 11) and the energy jump condition 
given by Eq. (3 13). The discontinuity of the velocity field, the velocity gradients, and 
the viscosity are treated by smoothing 


3.6 Smoothing of the Void Field 


The basic idea of all continuum surface tension (CST) models is the representation of 
surface tension as a continuous force per unit volume that acts in a neighborhood of the 
interface The non smoothed void field changes abruptly across the interface While 
calculating the surface tension force, applicable to a transitional region, first and second 
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order spatial derivatives of a are required These derivatives causes problem with CST 
formulation, if they are not smoothed out 


Here we calculate a smooth void field by convolving a discrete form of void fraction with 
a smooth kernel (a eighth-degree polynomial) This convolution results m a smooth void 
field, denoted by a, given by, 


a(x) = f a(x)T(x — x')dS (3 19) 

where r(z) is the kernel function and Clh is defined by a parameter e that essentially 
controls the smoothed interfacial thickness 


The choice of kernel is critical for the design of an accurate CST model The basic re- 
quirements for selecting a kernel function are as the following, 

(l) have compact support, 

(n) monotonically decreasing, 

(m) be radially-symmetric, 

(iv) be sufficiently smooth, i e , for some k > 3, kernel should be k times continuously 
differentiable, 

(v) have a normality property, i e , r(x, e) = 1, 

(vi) approach the Dirac delta function S(x) in the limit |D fc | — >■ 0 


In the present work we used a eighth degree polynomial kernel, denoted T, given by 


T(a:' - x) = { 


0 



\x — x | < e 
\x — x'\ > e 


(3 20) 


where A c is a normalization constant that ensures T(x,e) — 1 The kernel F(x, e) 
satisfies all the six properties listed above We have used e = 4 Sx, where 6x is the grid 
spacing The curvature, k, is calculated using the smoothed void field as, 


(3 21) 


k = —V n 
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where, 


Vet 
n = — — 
Va 


(3 22 ) 


3.7 Numerical Procedure 


The spatial discretization of governing equations is obtained using a staggered grid ar- 
rangement of Harlow and Welch [8] with scalars located at the cell centers and velocity- 
components located at the center of the cell faces The convection terms of the momentum 
equation are discretized by higher order methods using an ENO discretization scheme (see 
Chang et al [5]) The convection terms in energy equation are discretized by QUICK [14] 
The temporal discretization is described by a semi-implicit forward Euler method We 
begin a time cycle by solving the explicit energy equation m the bulk fluid phases 

9 n+1 = e n + il[_ v V(pc p 0) + V(/cW)] n (3 23) 

pc p 

After every time step with the help of new temperature field, thermal conductivity and 
specific heat at each cells are calculated The thermal properties and temperature field 
thus obtained are used to form the interfacial heat flux jump appearing m the mass source 
term The continuity and momentum equations are discretized m time as 



n dS + 


[ (- 

pi 


l) llq n+1 l 

Pg hlg 


n 


dS = 0 


(3 24) 


v n+l ^ v n _ W )n _ ^{ypn+1 + ( pg )n + y ^(y v + (W) T )] n + o(/cVd) n } (3 25) 

The velocity for the new time step is eliminated from these discrete equations and re- 
sulting pressure equation is solved by an iterative method based on a preconditioned 
conjugate gradient scheme of Van der Vorst [26] that has been implemented by Welch 
and Rachidi [28] Once the pressure at the new time level is obtained, the velocity at the 
new time level is found from the discrete momentum equations and the new density field 
is found from Eq (3 4) as 
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As discussed earlier, at this stage LVIRA algorithm is invoked The cells that are not mix- 
ture cells, or are not adjacent to mixture cells, do not require this treatment Once the 
new density field is found, the new void fraction may be calculated and the mixture cell 
mterfacial geometry can be updated 

The numerical scheme is based on the explicit time-advancement strategy The time 
step is determined from the limit imposed on a capillary wave traveling in an infinite 
medium [30] The wave is not allowed to travel more than half a cell width during a time 
step 



Chapter 4 


Results and Discussion 


The simulations were performed with variable and constant thermal properties of the 
fluid The variable properties are indicated in Table 4 1 The constant properties are 
evaluated at near-critical conditions of 373° C and 219 bar These properties are given 
m Table 4 2 The computational domain chosen for the present simulation is shown in 
Fig 3 1 Since the problem is symmetric, a computational domain with width A 0 /2 was 
chosen In the present simulation the numerical values of A 0 and A are 2 270 mm and 
0 2087 mm, respectively The height of the domain, H was set as A 0 with initial interface 
height, 5 set to 



We first present results of a convergence study undertaken to find an adequate grid res- 
olution For this, simulations were performed with three grid-meshes, viz , 90 x 180, 
180 x 360 and 360 x 720 Figure 4 1 shows the mterfacial shape at the instant of bubble 
release from the surface for the above-mentioned grid resolutions We observe that the 
180 x 360 and 360 x 720 produce nearly identical results, hence, we choose the 180 x 360 
grid for our subsequent simulations The same simulation run on 180 x 360 with the time 
step halved also gives nearly identical results as seen m Fig 4 2 Thus the 180 x 360 grid 
with the time step used in spatial convergence study provides results converged m time 
and m space, and it is used for all the subsequent simulations 

In the present simulation we considered water near the critical point and the vapor trans- 
portive thermal properties (thermal conductivity and specific heat of vapor) varying with 





Figure 4 2 Bubble interface for two different time steps Grid resolution is kept as 
180 x 360 
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temperature As mentioned earlier, the properties aie given in Table 4 1 The mterfacial 
properties and density are provided by Table 4 3 At the outset, we present numerical 
simulations showing growth of the bubble and its release pattern, with constant wall tem- 
perature boundary condition of 15°C superheat We used a data set of thermal properties 
(from Table 4 1) containing 10 values in the range of 15°C superheat and interpolated 
the properties to get intermediate values at any temperature corresponding to the cell 
temperature 


Temperature ° K 

k, W/mK 

c p , kJ/kgK 

p, Pas 

646 15 

0 5383 

3 520 x 10 2 

3 230 x 10" 5 

647 81 

0 2502 

0 4318 x 10' 2 

2 912 x 10~ s 

649 48 

0 2126 

0 2854 x 10* 

2 825 x 10~ 5 

651 15 

0 1918 

0 2227 x 10 2 

2 775 x 10" 5 

652 81 

0 1777 

0 1864 x 10 2 

2 743 x 10~ 5 

654 48 

0 1673 

0 1622 x 10 2 

2 720 x 10~ 6 

656 15 

0 1591 

0 1448 x 10 2 

2 704 x 10" 5 

657 82 

0 1524 

0 1316 x 10 2 

2 692 x 10" 5 

659 48 

0 1468 

0 1211 x 10 2 

2 683 x 10-* 

661 15 

0 1421 

0 1126 x W 

2 676 x 10" 5 


Table 4 1 Properties from NBS/NEC table for vapor phase. 


Substance 

k, W/mK 

c p , kJ/kgK 

p, Pas 

Liquid Water 

0 5454 

2 18 x 10 2 

4 67 x 10“ 5 

Water Vapor 

0 5383 

3 52 x 10 2 

3 23 x 10 -5 


Table 4 2 Properties used for constant property simulation 


a, N/m 

hi g , kJ/kg 

Qsat, K 

P saU P & 

p, kg/m 3 
Liquid water 

p, kg/m 3 
Water vapor 

0 07 x 10~ 3 

276 4 

646 15 

21 9 x 10 6 

402 4 

242 7 


Table 4 3 Interfacial properties and density 












(g) 
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Figure 4 3 shows the periodic bubble release patterns with growing interface and varying 
vapor volume The steady cyclic growth and release of vapor bubbles is usually called 
ebullition cycle Starting with a thin vapor layer, gradually the bubble grows bigger 
Finally, the bubble reaches a limit at which it detaches The detached bubble moves 
upwards and leaves the domain of interest The released bubble leaves a very thin vapor 
film near the nodal position on the wall on which it had rested prior to its departure 
Due to the vapor production at the interface the thin vapor, layer again grows in size 
The next bubbles tend to grow at the outer side walls, which are basically the anti-nodes 
according to the Taylor’s wave length 


The ebullition cycle or bubble release cycle is revealed through the sequence of varying 
interface boundary at different time instants. The mechanism for the repeating pattern 
can be described m the following way After the bubble is released from the surface, the 
vapor left behind experiences a downward force due to surface tension Subsequently this 
vapor packet is pushed down towards the film. The vapor packet impinges on the horizon- 
tal surface The surface tension induced flow promotes the movement of pressure gradient 
driven impinging vapor packet towards the symmetric side walls The vapor turns upward 
near the side walls to initiate an identical bubble release cycle Similar patterns of bubble 
release were observed by Son and Dhir [23] m their axisymmetric simulations of saturated 
film boiling 


Furthermore, the cyclic bubble release pattern can be discerned from Fig 4 4 and 4 5, 
which show the variation of space averaged wall heat flux with time for constant and 
variable thermal properties, respectively. For constant thermal properties (Fig 4 4), af- 
ter an initial transient, we observed a cyclic pattern with the time period of 0 27 sec 
For variable thermal properties (Fig 4 5), we found a smaller time period of 0 23 sec 
Figure 4 6 shows the variation of interface heat flux with time for both the cases, namely, 
constant and variable thermal properties It can be observed that with the variable ther- 
mal properties, a higher magnitude of energy is transferred from vapor to liquid This 
energy is responsible for conversion of liquid into vapor Figure 4 7 shows the variation of 
fractional vapor volume (vapor volume/ (vapor volume and liquid volume)) with time for 
both the cases Due to higher vapor generation, the bubble growth rate is faster for the 
case of variable thermal properties as compared with the constant thermal properties and 
thus, a smaller time period can be observed in the case of variable thermal properties. Son 
and Dhir [22] reported the bubble release time-period of the order of 0 1~0 15 seconds in 
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Results and Discussion 



Figure 4 4. Fluctuation of space averaged heat flux on the wall surface with constant fluid 

thermal properties 



Figure 4 5: Fluctuation of space averaged heat flux on the wall surface with variable fluid 
thermal properties 





Heat flux through liquid- vapor interface ( x 1 0 4 W/m 2 ) 
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Figure 4 6 Fluctuation of liquid-vapor interface heat flux with time 



Fractional vapor volume 
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Results and Discussion 



Figure 4 7- Variation of fractional vapor volume with time for constant and variable 
thermal properties 
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then numerical snnuinuoiis using water at 1 atm pressure The simulations by Son and 
Dim 1221 were at eomphshed with the consideration of large degree of superheats, of the 
ordei oi 100 ( Sint e m the piesent w'ork, the degree of superheat is 15° C, a relatively 
Luge tune peiiod ih expected As mentioned earlier, the vapor generation is higher for 
the \anable piopeitv simulation This observation is well supported by Fig 4 6 for heat 
flux thiough the liquid- vapor interface The interface area (m two-dimensional this is a 
cuned line) is used for this calculation However, the wall heat flux (Figs 4 4 and 4 5) for 
constant properties is higher then that for the variable properties There is an apparent 
contiadiction here This can be explained as the following Table 4 1 and 4 2 reveal that 
specific heat (c v ) of the constant property fluid is much higher than the average specific 
heat of the variable property fluid As the energy at any phase may be expressed m terms 
of the enthalpy (c p T) of that phase, the temperature of the vapor at constant property 
is much less than that of the vapor at variable property The higher temperature of the 
vapor (lor the case of variable property) sets up a higher temperature gradient at the 
interface entailing higher vapor generation. 


Figures 4 8(a) and (b) show the isotherms and mterfacial shape during the bubble growth 
at two different time instants Each isotherm almost follows the shape of the interface 
The spacing between the isotherms is m general uniform A deeper look mto Fig 4 8(a) 
reveals slightly nonuniform distribution of the isotherms in the primary bubble In the 
primary bubble, the isotherms are a little closely spaced near the neck due to narrowing 
of the cross-sectional area and the upward motion of vapor 


Figure 4 9 shows the temperature contours within the vapor film on the solid surface, near 
the location of minimum film thickness, for the degree of superheat 15° C The uniform 
spacing between the isotherms indicates that energy transfer from the wall, during film 
boiling, in this region, is governed primarily by conduction rather than convection This is 
further corroborated by Fig. 4 10, which shows the variation of peak heat flux value with 
the reciprocal of the minimum film thickness at various time instants The linear nature 
of variation indicates that heat flux is inversely proportional to vapor film thickness In 
addition, the nature of variation implies the dominance of conduction as mode of heat 
transfer near the location of the film with minimum thickness 


Figures 4 11 and 4 12 show the velocity vectors at different arbitrary time instants, 0 375 
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Figure 4 8 Isotherms for a superheat of A 6 = 15 °C in the vapor over the computational 
domain at time instants 0 375 s and 0 625 s respectively In (b) contour levels are spaced 
at a temperature difference of 1 °C for a range between 647 — 658°C 
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Figure 4.9 Isotherms for a superheat of A 6 = 15 °C in an enlarged region of Fig 11(a) 
near the location of minimum film thickness 
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Reciprocal of non dimensionaiized minimum film thickness /S) 


Figure 4 10 Variation of peak heat flux value with the reciprocal of minimum film thick- 
ness at various time instants 






eat of A 9 = 15° at time = 0 625 s 



Figure 4 13. Streamlines at time = 0 375 s 
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Figure 4 1 % \ ,v: itiun of heat flux on the wall surface at various time instants before 

bubble i* h on light symmetric boundary 


and 0 r,2f> snr respectively, for the case of variable thermal properties of the vapor m the 
computational domain Fur clarity, velocity vectors are plotted at the alternative grid 
points It can b.> seen that the upward movement of the interface is much stronger in the 
primal y bubble iegion (left side symmetric boundary m Fig 4 11) than m the region of 
secondary bubble (light side symmetric boundary in Fig 4 11) Figure 4 12 shows the 
velocity vectors just after the bubble has relased from right side symmetric boundary 
Furtheiinoie, the vapor injected into the primary bubble through the thm film induces 
vortex motion neni the interface which can be better observed in Figs. 4.13 and 4 14 
Such observations are quite similar to what have been found by Son and Dhir [22] 


The most important quantitative issue of film boiling is the wall heat transfer Although 
the wall temperature is constant m the present simulation, the heat flux from the wall 
varies over the wall surface. Figure 4 15 shows the distribution of wall heat flux for 15°C 
superheat at various time instants. We observe that the profile of the heat flux distnbu 
tion at a particular instant is just reciprocal to that of the vapor film thickness at that 
instant The heat flux from the surface is inversely proportional to the vapor film thick- 
ness at any wall location. Maximum heat transfer takes place at the location where the 
film is thinnest Under the bubble core, very little heat transfer takes place It is evident 
that the maximum and minimum heat flux corresponds to minimum and maximum vap 
film thickness, respectively. We observe the magnitude of the peak heat flux to increase 
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hum** 1 W ^ .uiH'i'ii. >'i hf.it tint on wall surface at time instants after the bubble release 

on n.‘h' 1 \ iii*" 1 '’ s nuuudan. 


nine. ( use, ,...*,„ ly. the peak shifts to the right. The figure signifies growth 

„j d (In light symmetric plane. This can be explained m the following way 

,\s ,i„. buhl,!,. „.** u* . the interface m the peak region moves upward, the, interface 
at ,„|,e, r,.,.,. ms tails to conserve the vapor volume. The vapor generation ,s unable to 
Otis.*, the ei.iiwit.il no. nt vapor volume by a large magnitude. Thus the film th f 
,s reduced over of the wall .surface except the region just underneath the bubble 

am* Furthermore, .he location at which the 61m is the thinnest, moves radially inward 
(considering .he pro,eclio„ of the spherical bubble . circular) on the sohd surface as the 
interface evoh es into a fully grown bubble The radially inward movement of the thrnnest 

layer explains, the shift of the peak 


Further, Fig. 4 10 shows the same trend ,n spatial heat Snx vanatmns “ " 

after the bubble release at the right symmetry boundary, wit pe ea flux 

to left symmetry boundary. We observe that there s is i st ^ ^ 

value just after the bubble release. The *■»«** temp eratnre gradient) at the 

interface can be defined as the nondimensiona 

wall m the following way 


Nu : 


*2 fAo/2 


r 

J 0 


A 


89 


($ w - 9 S at) dy 


dx 


(4 2) 


2/= 0 
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Figure 4 17 Nusselt number variation with time (with constant thermal properties) 


where A is the characteristic length Figures 4 17 and 4 18 show the evolution of the space 
averaged Nusselt number with time on the horizontal surface for constant and variable 
thermal properties, respectively A cyclic temporal variation (initial transients upto three 
cycles being ignored) of the space averaged Nusselt number has been observed We note 
that during the stage of bubble growth, the space averaged Nusselt number increases 
As the interface m the peak region moves upward, the vapor film at another location 
of wall becomes extremely thm and thereby entails enhanced heat transfer Just after 
the departure of the bubble, the heat transfer rate decreases leading to a sudden drop in 
Nusselt number which can be seen in the plot This is because the surface tension, acting 
as a restoring force, pulls down the interface at the location of the peak, and in turn, the 
interface m the valley region moves upward Thus the average film thickness over the ma- 
jor portion on the surface increases leading to a reduction m heat transfer Thereafter, as 
the next bubble ag ain starts growing, the Nusselt number increases and this cycle repeats 


Now we compare our time and space averaged Nusselt number values with the predictions 
made in earlier investigations The correlations due to Berenson [2] depict the Nusselt 
numbers constant with time and are not affected by the bubble dynamics unlike the 
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Figure 4.18 Nusselt number variation with time (with variable thermal properties) 


present numerical results. Berenson [2] has predicted the Nusselt number as 


f Gr Pr 1 
Nu = 0 42 < — r— ) 


1/4 


l 


T 


(4 3) 


Here Gr and Pr are the Grashof number and Prandtl number, respectively, defined by 


Gr 

Pr 


- Ll 


t%9>? ( Pi 


*9 

C PgPg 

hn 


— - 1 
P9 / 


(4 4) 


and 0 is the ratio of sensible heat to latent heat described as 


hfg 


(4 5) 


Klimenko [12] earned out a somewhat generalized analysis of film boiling on hori- 
zontal flat plates Employing a basic formulation, similar to that Berenson [2], Kilmenko 
developed a correlation that included data near critical pressure similar to the present 
numerical study According to his correlation, Nusselt number data for film boiling on an 
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npw.ml fat mg hon/ontal surface was expressed with m ±25 percent as 


\u - 1 fa) x Kr'Gi^Pr 1 ^, 

for Gr < 4 03 x 10 5 (4 6) 


\u = 2 16 x lO^Gr^Pr 1 / 3 ^, 


for Gr > 4 03 x 10 5 (4 7) 


where, 

fi = P 
= 0 895" 5 / 3 
f 2 = l 

= 0 7i p- 1 ' 3 


for 13 > 0 71 
for 0 < 0 71 
for 0 > 0 50 
for P < 0 50 


In Figs 4 17 and 4 18 we compare the span average Nusselt numbers caluculated 
with constant and variable thermal properties to the correlations of Berenson (Eq 4 3) 
and Klimenko (Eq 4 5 and 4.6) Equation 4 4 implicitly calculates the Nusselt number 
with the characteristic length and the thermal conductivity of the vapor evaluated at 
the wall tempeiature The correlations of Berenson and Klimenko both specify that the 
properties be taken at average film temperature ({T wa u+T sat )/2) The time-space average 
Nusselt numbers for both simulations and both correlations are listed m Table 4 4, from 
which it can be seen that there is reasonable agreement between all four values but the 
simulation with variable properties is closer to the correlations If the properties m the 
correlations are evaluated at the wall temperature, the Nusselt numbers for Berenson 
and Klimenko become 4 80 and 4 36 respectively This illustrates the shortcomings of 
attempts to account for the rapid variation m the fluid properties very near the critical 
point by evaluating them at a single temperature, whether m a correlation or a simulation 
The simulation with variable properties takes proper account of their nonlinear variation 
and defines the Nusselt number in a rational manner The discrepancies in the results 
between the values due to correlations and present analysis could also be attributed to 
the existence of wavelength shorter than the most dangerous wavelength as observed in 
the experiments of Hosier and Westwater [10] As also noted by Son and Dhir [22], this 
difference could be due to the inability of a two-dimensional model to take into account 
the additional variations pertaining to the interface position Such variations are expected 



Results and Discussion 


it) 


5 1 1 i 11 ' IVI< ^’ more efficient flow patterns for vapor removal and thus enhanced heat transfer 


' ( VtSi> 

Nu 

(computed) 

Nu (Berenson 
model) 

Nu (Klimenko 
correlation) 

% variation 
from Berenson 

% variation 
from Klimenko 

with \ a liable 
thermal 
; propei t ios 

3 88 

4 37 

4 20 

= - 11 2 

= - 76 

i with constant 
thermal 
i pi operties 

6 15 

4 37 

4 20 

= + 40 7 

= + 46 4 


Table 4 4 Comparison of predicted Nusselt number with that of Berenson and Klimenko 
con elation for constant and variable thermal properties 


Chapter 5 


Conclusion and Scope for Future 
Work 


Investigations of pool boiling have been carried out with the aim to employ numerical 
techniques to advance the understanding of bubble formation and the induced fluid mo- 
tions and consequently their influence on the instantaneous heat transfer characteristics 
from the wall surface This numerical investigation provided a nice insight into bubbly 
film boiling flows, yielding fascinating results The results of this work revealed 


• Successful predictions of bubble growth and heat transfer characteristics for film 
boiling on a horizontal surface using VOF method 


• The predictions have revealed unsteady periodic nature of bubble growth at the 
nodes and the anti-nodes 


• Film thickness and heat transfer coefficient are found to vary spatially and tem- 
porally during the growth of interface 



itude and location of maximum heat transfer from the wall surface are found 
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to he dependent on file 


magnitude and location of minimum film thickness 


* i’he Nusselt numbers for film boiling (with variable thermal properties) obtained 
m the present work aie only 11 2% and 7 6% lower than those obtained from Berenson’s 

model and Khmeims con elation respectively 


In the piesont work, fhe computations were performed under two-dimensional and sym- 
metric boundary conditions using a predefined separation between regions of bubble for- 
mation In future computations, the restrictive assumption of boundary conditions will 
lie relaxed find a fully three-dimensional prediction must be aimed for Future plans also 
include the use of our method m numerical studies of the transition region of the boiling 
curve, m which the problems associated with moving contact line will be addressed 
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Appendix A 


A.l List of Subroutines 

film.f - main piogram this calls all subroutines 

bc.f - boundary conditions 

delta.f - calculates dx,dy (cell size) 

density.f - calculates new density field (equation 3 26) 

energy.f - bulk phase energy equation (equation 3 23) 

property.f - After each iteration (in energy f) with the help of temperature value at 
the centei of the cell This subroutine calculates thermal conductivity and specific 
heat for each cell Here we used a data set of thermal properties (From NBS/NRC 
Table) containing 10 values at a range of 15 C superheat and interpolated it to get 
intermediate values at any temperature corresponding to cell center 

iface.f - interface reconstruction with the help of LVIRA algorthim 


indta.f - input data and also calculates normalization constant for kernel function 
init.f - initial conditions 

iters.f - bcgstab routine and supporting routines 

junpro.f - This subroutine computes explicit junction voids The junction voids are 
used for density and viscosity terms in the momentum equations 
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f f 1 ' Uc, itatos volume of liquid and vapor in an interface cell 

linM 'R f Pomts of the interface with the cell boundaries 

moved .f ■ ret. re of each cell after the calculation of new density Held 

out (it a f Minima data to data files 

pviter.f- bml.lr mnd,l,„l continuity equation that includes mass source terms and calls 

notation u .mint 

qmt.f - l his piogramme finds temperature of a 2-phase cell based on extrapolation and 
also dctei mines the mtei facial heat fluxes For calculation of interfacial heat flux 
terms (sourc v and sourcl), we requires exact value of thermal conductivity along the 
line joining fiom midpoint of linear interface of a cell (i,j) to the center of neighboring 

cell (lq.jq) qmt f calls a subroutine named point f for the calculation of exact value 

of t hernial conductivity 

point. f - 'I his subioutme finds out the point of intersection of the cell faces and the 
line th.it conects the center of the cell (iq,jq) with the mid point of phase interface 
Then, with the use of thermal resistance concept, equivalent thermal conductivity 
along this line is calculated, which is used m qmt f for the calculation of the source 

terms 

setup. f - initialize velocities and other cell variables and also sets the initial interface 

profile 

smooth.f - smoothing of void field with kernel function calculation (Eqn 3 19-3 22) 

solid. f - solution of energy equation m adjacent solid surface However, m the present 
simulation we set wall temperature in this subroutine 

I 

vexplt.f - momentum equations (Eqn 3 25) 
weight. f - weightings for smoothed gradient calculation 


A. 2 Flow Chart 
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